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Using molecular dynamics (MD) simulation, we investigate the mechanical response of silicon 
to high dose ion-irradiation. We employ a realistic and efficient model to directly simulate ion 
beam induced amorphization. Structural properties of the amorphized sample are compared with 
experimental data and results of other simulation studies. We find the behavior of the irradiated 
material is related to the rate at which it can relax. Depending upon the ability to deform, we 
observe either the generation of a high compressive stress and subsequent expansion of the mate- 
rial, or generation of tensile stress and densification. We note that statistical material properties, 
such as radial distribution functions are not sufficient to differentiate between different densities of 
amorphous samples. For any reasonable deformation rate, we observe an expansion of the target 
upon amorphization in agreement with experimental observations. This is in contrast to simulations 
of quenching which usually result in denser structures relative to crystalline Si. We conclude that 
although there is substantial agreement between experimental measurements and most simulation 
results, the amorphous structures being investigated may have fundamental differences; the differ- 
ence in density can be attributed to local defects within the amorphous network. Finally we show 
that annealing simulations of our amorphized samples can lead to a reduction of high energy local 
defects without a large scale rearrangement of the amorphous network. This supports the proposal 
that defects in amorphous silicon are analogous to those in crystalline silicon. 
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PACS numbers: 61.43.Dq, 61.72.Tt, 61.80.-x 



I. INTRODUCTION 

> 

0\ • Mechanical response of materials to ion irradiation has implications for many materials applications. Ion processing 
of silicon is of widespread fundamental and technological importance due to its central role in the semiconductor 
C*") . industry. The high dose, large atomic number wafer implants used today lead to a significant amount of amorphization, 
which must be removed by subsequent annealing. Ion implant induced stresses can lead to such problems as substrate 
bending, delamination and cracking, and anomalous diffusion of dopantsa. The creation and modification of stress in 
silicon by ion irradiation has been examined experimentally by several research groupsou. Experimental measurements 
have shown that Si expands by approximately 1.8% upon ion beam induced amorphizationd. Little atomic scale 
■ modeling of the amorphization process has been conducted, although many attempts have been made to simulate 
the structure of amorphous Si (a-Si). These simulations usually involve a rapid quench of liquid Si and suffer from 
I ■ drawbacks such as limited system size and short simulation times. The vast majority of these simulations result in 
a structure with similar properties to experimental a-Si, except for the fact that the density is typically around 4% 
greater than crystalline Si (c-Si)u. Here, we study the formation and structure of a-Si by direct MD simulation of the 
O ■ ion beam amorphization process, rather than by a simulation of quenching. 



II. MOLECULAR DYNAMICS SIMULATION SCHEME 



ft ■ 

We employ a minimal atomistic model of the radiation damage process and investigate stressing and deformation of 
the substrate. A full MD simulation incorporating a large enough section of the target material to completely contain 
the paths of many implanted ions, and to run for a realistic time is neither feasible, nor particularly desirable. What 
we wish to study is the response of a section of the material to the damage induced by the irradiation. In order to 
achieve this, we simulate a model system, that incorporates all the necessary features of the complete system. 

At the ion energies used, the material response is a bulk, rather than a surface effect, so it is not necessary to 
explicitly include the surface in simulations. Hence, periodic boundary conditions are applied in all directions to 
minimize finite-size effects. Energy dissipation and equilibration of pressure are controlled by coupling the system to 
an external bathu. This involves two parameters, tt and rp, that are time constants controlling the rate at which 
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temperature and pressure tend towards set values. This non-physical coupling gives us control over quenching and 
expansion of the system, but does limit the inferences we can draw about dynamic processes. We have analyzed the 
dependence of simulation results on the strength of these couplings in order to ensure that we do not introduce any 
systematic errors into our simulations. We assume that the ion dose is sufficiently low that the presence of ions in 
the material has a negligible effect and may be ignored. The effect of the irradiation is modeled by choosing primary 
knock-on atoms (PKAs) within the simulated region and assigning velocities to them from a distribution that can 
be dependent upon ion type, ion energy and the depth in the substrate. PKAs are produced by collisions with the 
implanted ion, but the majority are due to collisions of knock-on atoms in sub-cascades. Hence our PKA velocity 
distribution must take account of all knock-on atoms that are created by ion implantation. In order to account for 
the free surface, the material is allowed to expand or contract in the direction normal to the surface; expansion is 
not allowed in the in-plane directions. The passage of high energy ions through the simulated region are rare events, 
between which the material will relax as energy is dissipated, then enter an essentially constant state for a long period 
of time. As this constant state contributes nothing to the response of the material, we omit it and start another 
PKA at this time. By taking this approach, we can reduce the system size to the extent that a representative set of 
simulations can be conducted on a few workstations in a matter of days. 

III. SIMULATION DETAILS 

MD based ion implant simulations^ were used to to determine realistic PKA velocity distributions for various ion 
types and energies, at differing substrate depths; the resulting distributions were similar. The magnitude of PKA 
velocities could be approximately described by a \ 2 distribution. For an incident ion beam perpendicular to the 
surface, the in-plane distribution of PKA velocity vectors was uniform. The angle between the PKA velocity vector 
and the ion beam could be approximated by a Gaussian distribution with mean and standard deviation of 85° and 
10° respectively. 

PKA energies were chosen such that the size of displacement cascades was less than the system dimensions, to 
ensure that parts of the same damage cascade could not interact through the periodic boundaries. The selection 
of coupling constants for pressure and temperature control is critical; too long coupling times lead to excessive 
cpu requirements, whilst too short coupling times lead to unphysical system behavior, e.g., rapid quenching of hot 
disordered material will produce an excessive number of defects. Therefore, preliminary simulations were conducted 
with various combinations of PKA energy, system size and temperature and pressure coupling constants to investigate 
and minimize the sensitivity of the simulation on these parameters. 

The final simulations were conducted using PKA energies with a mean of 1 keV. All silicon atoms were coupled to 
a 300 K heatbath with a coupling constant, tt, of 1.0 ps; this was determined to be large enough that increasing it 
had no effect on simulation results. A new PKA was initialized as soon as the system temperature was quenched to 
below 320 K. Stress in the direction normal to the surface was relaxed to zero with a time constant, Tp, of 0.1, 1.0, 
or 10.0 ps. No single value of this parameter could be used as it effectively represents the ability of the material to 
locally deform, both elastically and through plastic flow. As such it is affected by the size of the ion path, ion mass 
and energy, the ion dose rate and the total ion dose. We therefore ran three sets of simulations spanning two orders of 
magnitude range in this parameter. This range is expected to contain most experimentally accessible conditions. In 
order to relate stress to change in atomic positions, a Young's modulus of 113.0 GPao was used for both crystalline and 
damaged Si; the value of this parameter is expected to vary during amorphization, but its exact value is not critical 
as variations will only affect the value of the the non-critical time constant, rpQ. The initial system was crystalline Si, 
with dimensions of approximately 50.0 A x 50.0 A x 50.0 A, containing 6240 atoms, at zero pressure. In order to be 
able to simulate a sufficiently large sample, we use a many-body empirical potential for silicorH To our knowledge, 
the derivation of the pressure virial for this potential has not been previously published. A method of obtaining the 
virial, which may be applied to any many-body potential, is given in Appendix |A|. 

IV. AMORPHIZATION OF SILICON DURING ION IRRADIATION 

Three sets of simulations are presented, the only difference between them being the magnitude of the pressure control 
coupling constant. In each case 9 separate simulations were conducted for each parameter set and the final results 
averaged. Data was recorded at the end of each quench, i.e., immediately prior to initializing a PKA, to generate 
radial distribution functions (RDFs), and distributions of bond angles, bond lengths, and atomic coordination. The 
simulation parameters chosen resulted in a mean time between PKA initialization of approximately 6.25 ps. 
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During irradiation, damage accumulates within the sample, resulting in the development of stress. As the simulation 
progresses, the material expands, or contracts to relieve the normal stress, as shown in Figs, [I] to ^. Anisotropic 
deformation occurs by flow of material to relieve the in-plane stress. The normal stress is relaxed to zero in the case of 
the two shorter Tp values and remains negative for the longer Tp. The stress would be expected to reach an equilibrium 
state for any value of Tp, given a sufficiently long simulation time, but cpu constraints prevented us from achieving 
this for the largest value of Tp. The in-plane stress is intially positive (compression) due to damage accumulation. 
Once a sufficient density of damage is present the sample is able to relax via flow of material and the stress becomes 
negative (tension). The in-plane negative stress is due to the outflow of material during the 'thermal spike' phase 
of radiation damage. As the sample cools during expansion a tensile state is reached, but material is not able to 
flow back tq-.relieve the stress. A similar response to ion irradiation has been previously observed in wafer curvature 
experimental. At long times, the damage appears to saturate and the material maintains a dynamic equilibrium, 
with a constant in-plane stress of around —0.3 GPa. The percentage expansion of the system is dependent upon the 
value of rp; smaller values (faster relaxation) result in expansion of the material, whilst the largest value gives an 
initial expansion followed by contraction. Due to cpu time constraints only the intermediate case was followed until 
the system reached an equilibrium condition. In this case .the material reaches a density around 1.9% less than c-Si, 



The distributions shown in Figs. [| to |7j were calculated for the final state of each amorphization simulation, i.e., they 
correspond to the final data points in Figs. [I] to [| The distribution of bond angles is shown in Fig. ||. To calculate 
the distribution, angles were weighted by the product of the Tersoff cut-off functions of the two bonds involved. The 
small peak at 60° corresponds to three-fold rings and five-fold coordinated atoms, the shoulder at around 80° is due 
to four fold rings, and the main peak is due to five and six fold rings. The distribution due to c-Si is also shown to 
illustrate the effect of damage and amorphization; the frequency is scaled by 1 /3 in order to allow the sharp peak to 
appear on the same plot as the a-Si results. The mean bond angle in the amorphized samples is around 2° less than 
the tetrahedral angle, with a standard deviation of approximately 19°. 

The distribution of bond lengths is shown in Fig. | . When interpreting this it is useful to recall that the Tersoff 
switching function acts between 2.7 and 3.0 A. The small peak between 2.9 and 3.0 A is therefore due to non- 
bonding (repulsive) interactions. Bond lengths become about 0.09 A greater than those in c-Si. The angle and length 
distributions due to the different relaxation rates are so similar that they cannot be easily distinguished. In order to 
count the coordination of atoms, we have to specify a criterion to decide if any pair of atoms are bonded. Based on the 
bond-length distribution we count all atoms with interaction distances within the mid point of the switching function 
(2.85 A) to be bonded. The resulting coordination distribution is shown in Fig. ^. Very few (~3%) undercoordinated 
atoms are produced during amorphization, but approximately 26% become five-fold coordinated, and occasional 
(~2%) six-fold coordination occurs. Even though the average coordination of atoms increases to approximately 4.27 
during amorphization, expansion is possible as bond-lengths also increase. 

Table [j] contains structural information on the samples produced by these simulations. The mean and standard 
deviation of the structural parameters are given, and the percentage expansion of the irradiated samples and potential 
energy per atom are given relative to c-Si at 300 K. The Tersoff potential gives a cohesive energy for Si of 4.63 eV/atom 
at zero K and 4.59 eV/atom at 300 K (this is in agreement with equi-partition of kinetic and potential energy, as 300 
K corresponds to a kinetic energy of 0.04 eV/atom). In the final damaged structures, the potential energy per atom 
is approximately 0.41 eV higher than that of c-Si at 300 K. This stored energy is approximately 8.9% of c-Si cohesive 
energy and is primarily due to bond angle distortions of around 19° as shown in Fig. || and Table |. 

The resulting radial distribution functions, shown in Fig. ^ are indistinguishable between the various simulation 
sets, but differ from the experimental data obtained by neutron diffraction experimentsEHj. The broader peaks in 
the simulated RDFs indicate that the structures contain a higher degree of disorder than the experimental samples. 
There is also a feature at around 3 A, where the artifact in the shape of the first minimum reflects the abrupt cutoff 
in the empirical potential. The first two peaks of the a-Si RDF correspond to those for c-Si showing the existence of 
short-range order and a tetrahedral like environment for each atom. There is little correlation between the positions 
of the other peaks, except for the third peak in the experimental data of Kugler et al.; this may be an indication of 
the presence of micro-crystallites in their sample. 



The damage simulations were carried out at room temperature with extremely limited times for structural relaxation. 
Therefore, the amorphized samples are expected to contain a large number of high energy defects; this is demonstrated 
by the difference between the experimental and simulated RDFs. In order to directly compare experimental and 
simulated samples, we have therefore carried out annealing simulations. The final structures from each set of damage 




V. ANNEALING 
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simulations were taken, and subjected to four separate anneals. Anneals were done at room temperature, between 
room temperature and the a-Si to c-Si transition temperature, between the transition temperature and the c-Si melting 
point, and at the melting point. The Tersoff potential over-predicts the c-Si melting temperature, giving a value of 
approximately 3000 KlI. Therefore we applied the following scaling between actual and simulated temperatures: 



Tersoff 



^Experiment ^Experiment — 300 K 

1.947 TExperiment ~ 283.994 K Tsxperiment > 300 K. 



This resulted in annealing temperatures of 300 K, 1121 K, 2194 K, and 3000 K, which approximately correspond 
to experimental temperatures of 300 K, 722 K, 1273 K, and 1687 K respectively. Anneals were continued until no 
further change in any of the measured structural parameters was observed. All simulation parameters were kept the 
same as those used for the damage simulations apart from the anneal temperature. After annealing, the samples were 
quenched to 300 K over a time of 40 ps before the final structural properties were measured. 

The structural changes during annealing are shown in the following plots for the rp=l ps damage simulation. The 
high temperature anneals resulted in large density changes as shown in Fig. [| After cooling to room temperature to 
remove the effect of thermal expansion, the final densities of the annealed samples corresponded to 2.1% and 1.9% 
expansion, and 0.4% and 1.0% contraction relative to c-Si respectively for the 300 K, 1121 K, 2194 K, and 3000 K 
anneals. The room temperature anneal gives very little structural modification in the time simulated; all distributions 
are indistinguishable from the unannealed sample. 

The 1121 K anneal results in modifications to bond angle, bond length, and coordination distributions and to 
the RDF, whilst the density is almost unchanged. This indicates that local defects are being removed without any 
global structural reorganization. This is in agreement with experimental studies which indicate that defects in a-Si are 
similar to defects in c-Si, i.e., they behave as interstitials and vacancies within the amorphous networkHj'til. Annealing 
therefore occurs through diffusion and annihilation of defects, rather than by a reorganization of the network. Fig. ^ 
shows the final distribution of bond angles; the mean increases by approximately 1°, to around 1° less than the c-Si 
angle, while the standard deviation of bond angles is reduced to 16°. Bond lengths are reduced by an average of 0.03 
A, to 0.06 A longer than those in c-Si, as shown in Fig. [l0| The anneal results in a reduction in 3-fold and 4-fold rings 
as evidenced by the reduction in the frequency of bond angles of 60° and 80° , and the reduction in 5-fold coordinated 
atoms. After annealing and quenching, approximately 1% of atoms are 3-fold coordinated, 83% of atoms are 4-fold 
coordinated, and 16% are 5-fold coordinated, as shown in Fig. |ll|. The calculated radial distribution function, shown 
in Fig. [L2L is in very good agreement with the experimental data, with the only discrepancy due to the short range 
of the potential cutoff. Annealing leads to an energy gain of approximately 0.13 eV per atom, to give an energy in 
the final relaxed structure 0.28 eV higher than that of c-Si. This stored energy is approximately 6.1% of c-Si cohesive 
energy and is primarily due to bond angle distortions of around 16° as shown in Fig. ^|. 

The 2194 K anneal causes partial melting of the amorphous material and subsequent reconstruction into a higher 
density structure. The final quenched structure contains more defects than the 1121 K annealed sample, but this may 
be in part due to the rapid quench from a semi-liquid state. We observe no indication of any recrystallization in the 
time simulated. The 3000 K anneal clearly involves complete melting of the sample, as the density increases to above 
that of crystalline Si. Again, the annealed sample contains more defects than the 1121 K annealed sample, but less 
than the un-annealed sample. 



VI. CONCLUSIONS 



We have investigated the structural response of crystalline silicon to ion-irradiation. We are able to simulate stress 
generation and structural relaxation with a relatively simple model. Simulation of continual radiation followed by 
annealing generates amorphous silicon with a low level of defects. Amorphous silicon prepared by various experimental 
methods usually has a density between 2% and 10% less than that of crystalline material, whereas atomistic simulations 
usually produce samples with a density around 4% higher than c-Si. The fact that a-Si prepared by simulations of 
quenching liquid Si has a higher density than c-Si is not surprising. Amorphous silicon has a structure similar to 
that of the liquid, which-is 5% denser than the crystal form at the melting pointQ. Therefore a defect free continuous 
random network (CRN)E9 structure for a-Si would be expected to result in a density greater than c-Si. 

We conclude that by conducting a direct simulation of radiation induced amorphization, we have produced the 
experimentally observed structure of ion beam amorphized Si. The structure differs from the proposed defect free 
CRN models, and differs from a-Si structures formed by simulations of quenching. Although this metastable structure 
is at a higher energy than CRN a-Si, it cannot be annealed to produce that structure, as transformation to c-Si, or 
melting will occur at a lower temperature. This leads us to question what is meant by 'amorphous' Si, as an 
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experimental sample termed amorphous may contain micro-crystals, CRN structures, and defects such as vacancies, 
interstitials, c-a boundaries, etc. In some sense the never experimentally observed high density structure can be 
regarded as true amorphous material, whereas other less dense materials must be regarded as containing defects at 
the very least, with a detailed structure that is preparation dependent. 
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APPENDIX A: DERIVATION OF THE VIRIAL FOR MANY-BODY POTENTIALS 



1. Introduction 



While the pressure virial is easy to calculate for a thermodynamic system described by pair potentials, the extension 
to systems modeled by many body potentials is not obvious. Here we derive the virial for Tersoff typeQ potentials, but 
note that the method is readily modified to other many body, or three body potentials so long as the configuration 
energy can be written down as a function of atomic coordinates. 

Following SmithE3 we write the pressure with an explicit dependence on volume. This is achieved by relating 
absolute atomic coordinates, r, to scaled coordinates, p through the volume, V, so that under isotropic expansion, or 
contraction of the system, the scaled coordinates of atoms are unchanged: 

r = V 1/3 p. (Al) 

The pressure, P, is then given by: 

P = NkT/V - $/3V, ( A2) 

where the virial, is: 

$ = W(dU([V 1/3 p] N )/dV). (A3) 



2. Application to the Tersoff Potential 



The configuration energy, U , is given by: 

N 

^EEmO, ( A4 ) 

i— 1 j<i 



where Uij is the energy of the bond between atoms i and j, as a function of r , the set of atomic coordinates. The 
functional form is: 



Uij = f(rij)[V R (rij) - hjVAinj)], 



(A5) 



where fifij) is a cutoff function that makes the interactions short ranged, Vn{rij) and VA{Tij) are pair terms, r\j is 
the bond length, and bij is a many body term. As bij is a function of the positions of neighboring atoms of i and j, 
these atoms experience a force due to the i-j bond. Let t he set of atom s in volved in the calculation of tty, i.e. atoms 
i, j and neighbors, be denoted by Sij. Substituting Eq. flA^ ) into Eq. ( |A3| ) gives: 



A' 



$ = 3V ^(EE^(^ 1/ V] Ar )/w). 



(A6) 



=1 j<i 



Since from Eq. (Al) 
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dr = ly- 2/3 pdV, 



this can be rewritten as: 



A" 



i—1 j<i 



(A7) 



(A8) 



Explicitly writing the dependence on interacting atoms gives 

N 



1=1 j<i fceSij 



(A9) 



or, in terms of atomic positions: 



iY 



* = <£EE -#■'*> 

i=l j'<i keSij 



(A10) 



where is the force on atom k due to the bond between atoms i and j. The atomic positions, r^, are local to the bond 
being considered, e.g., periodic boundary conditions are accounted for before the potential and force calculations. All 
local atomic coordinates can in fact be translated to be relative to the position of atom i for each bond, so the force 
on atom % can be neglected: 



N 



* = <EEE -#■'*>• 

i=l j<i fees',. 



(All) 



If the interaction potential is pairwise, i.e., bij depends only on atoms i and j, Eq. (All) reduces to the pair virial: 



JV 



*=1 j<i 



Individual elements of the pressure tensor, can be obtained by replacing Eq. (Al) with: 

r = Lp, 

where L is a 3 x 3 matrix with determinant V. For example, by using: 














1 











1 



(A12) 



(A13) 



(A14) 



one obtains: 



A' 



••EE E 

i=i j<i kes' . 



(A15) 
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FIG. 1. Variation in normal stress during amorphization. 
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FIG. 2. Variation in in-plane stress during amorphization. 
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FIG. 3. Change in height of simulation cell during amorphization. 
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FIG. 9. Final bond angle distributions in annealed samples. 
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FIG. 10. Final bond length distributions in annealed samples. 
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FIG. 11. Final distributions of atomic coordination in annealed samples. 
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12. Final radial distribution functions in annealed samples. 



TABLE I. Structural properties of simulated samples. 



Sample 



Bond Length (A) Bond Angle (° 



Coordination 



Expansion (%) 



Energy (eV/atom) 



c-Si 

Amorphized: 

TP =0.1 pS 
Tp — 1.0 PS 

t p =10.0 ps 
Annealed: 
300 K anneal 
1121 K anneal 
2194 K anneal 
3000 K anneal 



2.36±0.05 

2.44±0.13 
2.45±0.14 
2.45±0.14 

2.44±0.13 
2.42±0.11 
2.43±0.12 
2.43±0.13 



109.42±2.77 

107.69±19.20 
107.58±19.44 
107.44±19.76 

107.60±19.40 
108.32±16.16 
107.95±17.60 
107.77±18.26 



4.0 

4.22±0.55 
4.27±0.54 
4.32±0.54 

4.26±0.54 
4.15±0.40 
4.23±0.45 
4.26±0.48 



0.0 

4.47 
1.91 
-0.28 

2.12 
1.91 
-0.36 
-1.04 



0.0 

0.41 
0.41 
0.38 

0.37 
0.28 
0.30 
0.31 



11 



